metric_isomap <- metricML(x, s = s, k = k, radius = radius, method = method, annmethod = annmethod, eps = 0, distance = distance, treetype = treetype, searchtype = searchtype)
summary(metric_isomap)

# Comparison of Isomap embedding plot
par(mfrow = c(1, 2))
plot(metric_isomap$embedding, col = viridis::viridis(24)) # metricML
# plot(e@data@data, col = viridis::viridis(24)) # embedding from dimRed, same as metricML()
emb_isomap <-
  feather::read_feather(here::here("data/embedding_isomap_1id336tow.feather"))
plot(emb_isomap$`0`, emb_isomap$`1`, col = viridis::viridis(24)) # embedding from megaman, radius = 0.4

# Use adjacency matrix as input for metricML and dimRed, then the embedding should be close, as well as the Riemannian metric
metric_isomap$adj_matrix[1:10, 1:10]# renormalized adjacency matrix

metric_isomap$weighted_graph %>% is.connected()

riem_isomap <- metric_isomap$rmetric
riem_isomap[,,1] %>% 
  # isSymmetric()
  matrixcalc::is.positive.definite() # FALSE

np = import("numpy")
H_isomap <- np$load(here::here("data/hmatrix_isomap_1id336tow.npy"))
H_isomap[1,,, drop=TRUE] %>% 
  matrixcalc::is.positive.definite() # TRUE

pd <- rep(NA, N)
for (i in 1:N) {
  pd[i] <- 
    riem_isomap[,,i] %>%  # R
    # H_isomap[i,,, drop=TRUE] %>% # Python
      # isSymmetric()
      matrixcalc::is.positive.definite() # FALSE
}
pd %>% sum

# determinant(riem_isomap[,,i])
# mat <- riem_isomap[,,i]
# eigen(mat)

# Function for graph Laplacian
# Input: W: N*N weight matrix for the neighborhood graph, radius: bandwidth parameter
# Output: L: N*N graph Laplacian matrix
Laplacian <- function(W, radius, lambda = 1){
  
  D <- Matrix::Diagonal(x = rowSums(W)^(-lambda)) # inverse of a diagonal matrix
  W1 <- D %*% W %*% D
  D1 <- Matrix::Diagonal(x = 1 / rowSums(W1)) # inverse of Tn1
  L <- 1 / (radius^2 / 4) * (D1 %*% W1 - Matrix::Diagonal(nrow(W))) # c=1/4 for heat kernel, depending on the choice of weights, P15 footnote
  
  return(L)
}

# Function for Riemannian metric for each point
# The Riemannian metric and its dual associated with an embedding Y. 
# The intrinsic dimension d
# The Riemannian metric is currently denoted by G, its dual metric by H, and the Laplacian by L. 
# G at each point is the matrix inverse of H.
riemann_metric <- function(Y, laplacian, d){
  
  # TODO: add dimension check for all inputs
  
  H <- array(NA, dim = c(d, d, nrow(Y)))
  
  for (i in 1:d) {
    for (j in i:d) {
      yij <- Y[, i] * Y[, j]
      H[i, j, ] <- as.vector(0.5 * (laplacian %*% yij - Y[, i] * (laplacian %*% Y[, j]) - Y[, j] * (laplacian %*% Y[, i])))
      H[j, i, ] <- H[i, j, ] # symmetric matrix
    }
  }
  
  ## Pseudo inverse of H gives the final embedding metric h
  ## Array H corresponds to \tilde{H} in Step 4(a)
  ## The embedding metric H is the pseudo inverse of \tilde{H}
  for (i in 1:nrow(Y)) {
    Hsvals <- eigen(H[ , ,i])
    Huu <- Hsvals$vectors
    Hvv <- Hsvals$values[1:d] # top d largest eigenvalues, already sorted in decreasing order
    Hvv1 <- diag(x = 1 / Hvv)
    H[ , ,i] <- Huu %*% Hvv1 %*% t(Huu)
    H[, , i] <- 0.5 * (H[, , i] + t(H[, , i])) # fix not symmetric issue
  }

  return(H)
}

embed_den <- as_tibble(cbind(x = x[,1], y = x[,2], z = f)) %>%
  drop_na()
# library(akima)
pixel <- 100
lenbreak <- 5
akima.spl <- interp(embed_den$x, embed_den$y, embed_den$z, nx=pixel, ny=pixel, linear=FALSE)
p.zmin <- min(embed_den$z,na.rm=TRUE)
p.zmax <- max(embed_den$z,na.rm=TRUE)
breaks <- pretty(c(p.zmin,p.zmax), lenbreak)
# hdrcde
hdrscatterplot(embed_den$x, embed_den$y, noutliers = 10)

hdrinfo <- hdr.2d(embed_den$x, embed_den$y, prob = c(50, 95, 99))
plot.hdr2d(hdrinfo)
# Contour
contour(akima.spl, main = "smooth interp(*, linear = FALSE)", col = viridis(length(breaks)-1), levels=breaks, add=TRUE)
points(embed_den, pch = 20)

filled.contour(akima.spl, color.palette = viridis,
               plot.axes = { axis(1); axis(2);
                 title("smooth interp(*, linear = FALSE)");
                 points(embed_den, pch = 3, col= hcl(c=20, l = 10))})